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ABSTRACT 

We  consider  the  solution  of  the  linear  systems  of  algebraic  equations  which 
arise  from  elliptic  finite  element  problems  defined  on  composite  meshes.  Such 
problems  can  systematically  be  built  up  by  introducing  a  basic  finite  element 
approximation  on  the  entire  region  and  then  repeatedly  selecting  subregions,  and 
subregions  of  subregions,  where  the  finite  element  model  is  further  refined  in 
order  to  gain  higher  accuracy.  We  consider  conjugate  gradient  algorithms,  and 
other  acceleration  procedures,  where,  in  each  iteration,  problems  representing 
finite  element  models  on  the  original  region  and  the  subregions  prior  to  further 
refinement  are  solved.  We  can  therefore  use  solvers  for  problems  with  uniform  or 
relatively  uniform  mesh  sizes,  while  the  composite  mesh  can  be  strongly  graded. 

In  this  contribution  to  the  theory,  we  report  on  new  results  recently  obtained 
in  joint  work  with  Maksymilian  Dryja.  We  use  a  basic  mathematical  frame  work 
recently  introduced  in  a  study  of  a  variant  of  Schwarz'  alternating  algorithm.  We 
establish  that  several  fast  methods  can  be  devised  which  are  optimal  in  the  sense 
that  the  number  of  iterations  required  to  reach  a  certain  tolerance  is  independent 
of  the  mesh  size  as  well  as  the  number  of  refinement  levels.  This  work  is  also  tech- 
nically quite  closely  related  to  previous  work  on  iterative  substructuring  methods, 
which  are  domain  decomposition  algorithms  using  non-overlapping  subregions. 

1.  Introduction.  In  this  paper,  we  consider  the  solution  of  the  large  linear  systems 
of  algebraic  equations  which  arise  when  working  with  elliptic  finite  element  approximations 
on  composite  meshes.  In  this  contribution  to  the  theory,  we  report  on  new  results  recently 
obtained  in  joint  work  with  Maksymilian  Dryja. 

Finite  element  models  on  composite  meshes  can  systematically  be  built  up  inside  a 
frame  work  of  conforming  finite  elements;  cf.  Ciarlet  [3].  We  do  so  to  be  able  to  use 
a  number  of  technical  tools  which  are  available  primarily  in  the  conforming  case.  We 
begin  by  introducing  a  basic  finite  element  approxinmation  on  the  entire  region  and  we  then 
repeatedly  select  subregions,  and  subregions  of  subregions,  where  the  finite  element  model 
is  further  refined.  We  solve  the  resulting  linear  system  by  using  an  iterative  method  such  as 
the  conjugate  gradient,  Richardson  or  Chebyshev  method.  We  accelerate  the  convergence 
by  solving  so  called  standard  problems.  These  correspond  to  the  finite  element  models  on 
the  original  region,  prior  to  any  refinement,  and  those  on  the  subregions  prior  to  further 
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refiiiement.  We  can  thus  use  solvers  for  problems  with  uniform  or  relatively  uniform  mesh 
sizes,  while  the  composite  mesh  can  be  strongly  graded. 

This  approach  offers  a  number  of  advantages.  An  existing  code  can  be  upgraded,  in 
order  to  increase  the  accuracy  locally,  without  a  radical  redesign  of  the  data  structures  etc., 
since  we  can  use  the  old  code  to  solve  one  or  several  of  the  standard  problems.  Issues  of  data 
structures  and  geometry  are  generally  simpler  if  we  design  programs  for  composite  mesh 
problems  in  terms  of  simpler  standard  problems.  The  use  of  simple  standard  problems  also 
tends  to  improve  the  performance  of  the  programs  on  vector  machines.  Finally  we  note 
that  if  each  of  the  standard  problems  has  appreciatively  fewer  degrees  of  freedom  than  the 
composite  model,  then  we  might  benefit  from  solving  a  nimiber  of  smaller  problems  rather 
than  one  large  one.  In  other  words,  we  might  view  our  approach  as  a  divide-and-conquer 
strategy. 

We  study  two  families  of  algorithms  which  we  call  multiplicative  and  additive  respec- 
tively. The  so  called  additive  algorithms  are  partictdarly  well  suited  for  parallel  computing 
in  that  in  each  iteration  step  we  can  simultaneously  solve  all  the  standard  problems  on  the 
the  different  levels  of  refinement.  Synchronization  between  the  processors  is  required  once 
in  each  iteration,  namely  when  the  residual  is  computed  and  assembled.  One  or  a  group 
of  processors  can  therefore  be  assigned  in  a  straightforward  way  to  each  of  the  standard 
problems.  As  we  will  see.  we  can  view  the  multiplicative  algorithms  as  preconditioned  con- 
jugate gradient  methods,  while  the  additive  variant  involves  the  solution  of  a  transformed 
equation  with  the  same  solution  as  the  original  one  using  no  further  preconditioning. 

The  systematic  study  of  the  methods  under  consideration  goes  back  at  least  to  1983, 
when  the  Fast  Adaptive  Composite  (FAC  )  method  was  introduced  by  McCormick  [12].  Is- 
sues related  to  the  implementation  of  this  method  on  parallel  computers  led  to  the  introduc- 
tion of  Asynchronous  FAC  (AFAC)  methods  [8]  a  few  years  later.  Nimierical  experiments 
have  now  been  reported  for  model  cases  and  recently  Richard  Ewing,  Steve  McCormick  and 
others  have  begun  to  test  the  algorithms  for  more  difficult  problems  arising  in  industry. 
The  convergence  of  the  FAC  method  is  discussed  in  McCormick  et  al.  [12],  [13]  under  a 
certain  additional  regularity  assumption.  An  important  contribution  to  the  theory  is  given 
by  Bramble,  Ewing,  Pasciak  and  Schatz  [2],  who  outlined  a  proof  of  optimality  for  a  two 
level  FAC  algorithm.  They  require  no  additional  regularity  in  their  proof.  In  recent  papers 
Bj0rstad  [1]  and  Mandel  and  McCormick  [11],  develop  a  theory  for  both  multiplicative  and 
additive  algorithms,  using  two  levels.  In  particular,  they  obtain  interesting  results  concern- 
ing the  relationship  between  the  spectra  of  the  two  methods.  In  a  recent  paper.  Mandel 
and  McCormick  [10]  establish  optimality  for  a  multi-level  AFAC  algorithm  for  a  special 
model  problem. 

Our  study  begim  with  the  discovery  that  FAC  and  AFAC  methods  have  a  structure 
quite  similar  to  that  of  the  classical  Schwarz  procedure,  see  Schwarz  [14],  and  an  additive 
variant  thereof  recently  considered  in  Dryja  [4]  and  Dryja  and  Widlund  [5].  Our  earlier 
work  was  in  turn  inspired  by  a  recent  paper  by  P.-L.  Lions  [9]  in  which  a  variational 
frame  work  for  the  classical,  miiltiplicative  Schwarz'  method  is  developed  for  continuous 
elliptic  problems.  We  were  able  to  establish  a  rate  of  convergence  which  is  independent  of 
the  number  of  degrees  of  freedom  as  weU  as  the  number  of  subregions  for  a  special  kind 
of  additive  Schwarz  algorithm,  see  Dryja  [4]  and  Dryja  and  Widlund  [5].  In  our  view, 
the  central  theoretical  issue  of  the  iterative  refinement  algorithms,  which  are  discussed  in 
this  paper,  is  similar  to  that  of  Schwarz  type  algorithms  namely  the  design  and  study  of 
algorithms  for  which  the  rate  of  convergence  is  independent  of  the  number  of  subproblems, 
i.e.  the  ntmiber  of  refinement  levels  as  weU  as  and  the  mesh  sizes. 

In  section  2,  we  introduce  the  finite  element  problems  on  composite  meshes,  certain 


projections  and  our  algorithms  in  their  basic  form. 

In  Section  3,  we  consider  a  multiphcative  algorithm  for  the  multi  level  case  and  show 
that  its  rate  of  convergence  is  optimal. 

In  Section  4,  we  introduce  several  multilevel  additive  algorithms  and  provide  a  number 
of  bounds.  The  principal  result  is  that  one  of  these  methods  has  a  rate  of  convergence  which 
is  independent  of  the  ntunber  of  refinement  levels,  as  well  as  the  mesh  sizes. 

2.    Composite  Finite  Element  Problems  and  Basic  Iterative  Methods.    We 

consider  linear,  self  adjoint,  elliptic  problems  discretized  by  finite  element  methods  on  a 
bounded  Lipschitz  region  fl  in  i?".  To  simplify  the  presentation,  we  assume  that  the  differ- 
ential operator  is  the  Laplacian  and  that  we  use  continuous,  piecewise  linear  finite  elements. 
However,  almost  all  of  our  results  can  be  extended  immediately  to  general  conforming  finite 
element  approximations  of  any  self  adjoint  elliptic  problem  which  can  be  formulated  as  a 
minimization  problem.  The  continuous  and  discrete  problems  are  of  the  form 

a{u,v)  =  f{v),  'i  V  e  V, 

and 

a{uh,Vh)  =  f{vh),  y  Vh  e  V'',  (1) 

respectively.  The  spaces  V  and  V^  are  defined  in  the  next  few  paragraphs.  The  bilinear 
form  is  defined  by 

a{u,  v)  =        Vti  •  Vv  dx  . 
Ju 

This  form  defines  a  semi-norm  \u\h^  =  (a(u,u))^/"  in  H^{Q,). 

To  simpHfy  the  presentation,  we  assume  that  we  have  a  zero  Dirichlet  boundary  con- 
dition on  dn,  the  boundary  of  Q.  The  space  V  is  thus  Hq{Q).  We  note  that  we  equally 
well  could  have  considered  an  inhomogeneous  Dirichlet  problem.  The  space  V^  is  defined 
on  a  composite  triangulation,  which  is  possibly  the  result  of  a  large  number  of  successive 
refinements.  The  triangulation  of  17  is  given  in  the  following  way. 

We  first  introduce  a  relatively  coarse  triangidation  of  Q,  also  denoted  by  Qi,  and 
denote  the  corresponding  space  of  finite  element  functions  by  V'*' .  We  can  thinlc  of  this 
space  as  having  a  relatively  uniform  (or  uniform)  mesh  size  hi.  Let  fio  be  an  area  where  we 
wish  to  increase  the  resolution.  We  do  so  by  subdividing  the  elements  and  introducing  an 
additional  finite  element  space  V^^.  We  assure  that  the  resulting  composite  space  V^^  + 
V'^^  is  conforming  by  having  the  functions  of  V''-  vanish  on  d^2-  We  repeat  this  process  by 
selecting  a  subregion  ^3  of  ^2  and  introducing  a  further  refinement  of  the  mesh  and  finite 
element  space  etc..  We  denote  the  the  resulting  nested  subregions  and  subspaces  by  fl,  and 
V^'  respectively.  Throughout,  we  have  Qi  C  fii-i  and  V'''-'  (1  H^(Qi)  C  V^'  r  H^iQi), 
i  =  2,. .  .,k.  The  composite  finite  element  space  on  the  repeatedly  refined  mesh,  is 

yh  _  yhi  j^  yh2  ^  _  _  _  ^.  yf^i 

We  assume  that  all  the  elements  are  shape  regular  in  the  sense  that  there  is  a  uniform 
bound  on  hfc/pK-  Here  hf{  and  pf^  are  the  diameter  and  the  radius  of  the  largest  inscribed 
sphere  of  any  element  K,  respectively.  Our  bounds  in  the  theory  developed  below  also 
depend  on  the  shape  of  the  subregions  Qi.  Thus  in  order  for  our  proofs  to  work,  we  cannot 
allow  the  sets  f2,_i  \Qi  to  become  arbitrarily  thin  in  comparison  with  the  diameter  of  r2,_i. 
We  also  assume  that  the  area  of  any  triangle  on  level  i  can  be  bounded  by  const.  q-''\  j  <  i 
times  the  area  of  the  triangle  on  level  j  of  which  it  is  a  part.  Here  q  is  a.  constant  <  1. 


The  finite  element  problem  is  defined  by  equation  (1)  and  the  corresponding  stiffness 
miatrix  can  conveniently  be  computed  by  using  a  process  of  subassembly.  Introducing 
subscripts  to  indicate  the  domain  of  integration,  we  write 

a{u,  v)  =  au^\Q^{u,  v)  +  aQ^\Q^{u,  v)  +  . . .  +  an^(u,v)  . 

The  stiffness  matrices  corresponding  to  the  regions  fli  \  ri,+i,  i  <  k  -  1,  and  Clf.  are 
computed  by  working  with  basis  functions  related  to  the  mesh  size  hi  .  The  quadratic 
form  corresponding  to  the  composite  stiffness  matrix  is  the  simi  of  the  quadratic  forms 
corresponding  to  Qj-  \  fii+i  and  Clf..  In  our  algorithms,  we  use  solvers  for  the  same  elliptic 
problem  on  the  subregions  Q,i  ,  i  =  1, ...,  k,  and  the  relatively  uniform  meshes  corresponding 
to  hi  .  The  corresponding  stiffness  matrices  can  in  fact  be  obtained  at  a  small  extra  cost 
during  the  assembly  of  the  composite  mesh  model.  When  we  refine  a  finite  element  model 
locally,  the  modified  stiffness  matrix  is  obtained  by  replacing  the  quadratic  form  associated 
with  the  subregion  in  question  by  the  one  corresponding  to  the  refined  model  on  the  same 
subregion.  It  is  therefore  relatively  easy  to  design  a  method  which  systematically  generates 
the  stiffness  matrices  for  all  the  standard  problems  necessary  while  the  stiffness  matrix  of 
the  composite  model  is  computed. 

The  fundamental  building  blocks  of  our  algorithms  are  the  Pj,  i  <  _/',  the  projections 
onto  the  spaces  V^'  f]  Hq(Q.j).  We  note  that  if  j  >  i,  then  we  solve  a  problem  on  17^  with 
a  coarser  mesh  than  if  V'^'  were  used.  The  projection  Pj,  i  <  j,  is  defined,  in  terms  of  the 
unique  element  of  V^'  D  Hq(Q.j),  which  satisfies 

aiP'-Vh,<Ph)  =  a{vh,<Ph),  ^<l>h  e  V'''.  (2) 

We  now  introduce  the  multiplicative  and  additive  algorithms.  Since  no  further  effort 
is  involved,  we  develop  a  framework  which  is  also  usefiil  in  other  contexts  such  as  the 
study  of  algorithms  of  Schwarz  type;  cf.  Dryja  [4],  Dryja  and  Widlund  [5]  and  Lions 
[9].  We  begin  by  considering  the  case  of  two  subspaces  Vi  and  V2  and  the  multiplicative 
(sequential)  algorithms.  In  a  first  fractional  step,  we  find  a  correction  5iu"  G  V'l  of  the 
current  approximation  u"  by  solving 

a(<5iu",u)=  /(t;)-a(u",u)  =  a{u'  -  u",  u),  V  t;  G  Fi. 

Here  u*    e   Vi  +  V2  is  the  solution  of  the  given  problem.     The  calculation  of  a  second 
correction  52u"  e  V2  completes  the  (n  +  l)th  step. 

0(^2"",  w)  =  fiv)  -  a(u"  +  ^iu",i;)  =  a{u*  -  (u"  +  *iu"),u),  Vu  G  V;  . 

As  shown  in  Lions  [9],  it  is  easy  to  see  that 

5iu"  =  Pi(u'  -u") 

<52u"  =  P2(w*-u"  -^iu")  =  P.{I  -P){u*  -n") 

and  thus  the  error  propagates  as 

Here  P;  and  P2  are  the  orthogonal  projections  associated  with  the  bilinear  form  a(-,  •)  and 
Vi  and  V2,  respectively. 


We  can  thus  view  this  algorithm  as  a  simple  iterative  method  for  solving 

(Pi  +  P2  -  P2Pi)Uh  =  9h, 

with  an  appropriate  right  hand  side  gn-  We  note  that  this  operator  is  a  polynomial  of  degree 
two  and  thus  not  ideal  for  parallel  computing  ,  since  two  sequential  steps  are  involved  . 
If  we  use  more  than  two  subspaces  and  therefore  more  projections  this  effect  is  further 
pronounced.  The  basic  idea  behind  the  additive  form  of  the  algorithm  is  to  work  with 
a  simplest  possible  polynomial  in  the  projections.  With  k  subspaces  we  thus  solve  the 
equation 

PuH  =  {Pl+P2+...  +  Pk)Uh  =  g'h,  (3) 

by  an  iterative  method.  K  we  can  show  that  the  operator  P  is  symmetric  and  positive 
definite,  the  iterative  method  of  choice  is  the  conjugate  gradient  method  at  least  on  a 
computer  with  conventional  architecture.  We  must  also  make  sure  that  this  equation  has 
the  same  solution  as  equation  (1),  i.e.  we  must  find  the  correct  right  hand  side.  Since  by 
equation  (1),  we  have 

a{uh,<f>h)  =  f{<t>h)  , 

we  can  construct  the  right-hand  side  gj,  by  solving  this  equation  restricted  to  aU  the  different 
subspaces  and  adding  the  results.  It  is  similarly  possible  to  apply  the  operator  P  of  equation 
(3)  to  any  given  element  of  V^  by  applying  each  projection  Pi  once.  Most  of  the  work  ,  in 
particular  that  which  involves  the  individual  projections,  can  be  carried  out  in  parallel. 

It  is  well  known  that  the  number  of  steps  required  to  decrease  an  appropriate  norm  of 
the  error  of  a  conjugate  gradient  iteration  by  a  fixed  factor  is  proportional  to  ^/K,  where 
K  is  the  condition  nimiber  of  P;  see  e.g.  Golub  and  Van  Loan  [7].  We  therefore  need  to 
establish  that  the  operator  P  of  equation  (3)  is  not  only  invertible  but  that  satisfactory 
upper  and  lower  bounds  on  its  eigenvalues  can  be  obtained. 

We  end  this  section  by  discussing  a  symmetric  version  of  the  multipUcative  algorithm. 
The  projections  generally  do  not  commute  and  if  we  wish  to  have  a  symmetric  expression  in 
the  operators  P;,  which  allows  us  to  accelerate  the  convergence  by  the  standard  conjugate 
gradient  method,  we  have  to  use  additional  fractional  steps.  In  the  case  of  two  subspaces. 
we  can  solve  the  first  problem  again.  Since  (/  -  P2)^  =  (7  -  P2)  we  can  write  the  restdting 
operator  as 

I  -{I  -  Pi)(7  -  P2)(/  -  Pi)  =  Pi  +  P2  -  P1P2  -  P2P1  +  PiP:Pi  =  /  -  TjT;, 

where  T2  =  (7  -  Pi  )(7  -  P2 ).  The  corresponding  operator  for  k  subspaces  involves  '2k  -  1 
fractional  steps  and  has  the  form 

7  -  nn,  where  T,  =  (7  -  Pi)(7  -  P2)  •  ■  -(7  -  Pfc). 

The  error  propagation  operator  for  the  basic  multiplicative  algorithm  is  Tj^.  while  the 
convergence  rate  of  the  symmetrized  algorithm  can  be  bounded  in  terms  of  the  condition 
number  of  7  —  T}.T^.  We  note  that  a  bound  on  the  spectral  radius  of  T/.  is  obtained 
immediately  from  the  spectral  bounds  on  the  symmetric  operator. 

3.  The  Multiplicative  Algorithm.  In  this  section,  we  establish  the  following  result. 
We  note  that  our  analysis  resembles  that  of  Bramble  et  al.  [2)  in  the  case  of  k  =  2. 

Theorem  1.  The  symmetrized,  m.ultiplicative  iterative  refinement  algorithm  based  on 
the  projections  P;,  J  <  k,  has  a  condition  number  which  is  independent  of  k  arid  the  number 


of  degrees  of  freedom.  The  spectral  radius  of  the  basic  multiplicative  algorithm  is  bounded 
by  a  constant  which  is  uniformly  less  than  one. 

We  begin  by  designing  a  preconditioner  b{u,v)  for  the  finite  element  problem  on  the 
space 

yh  _   yh,    _|_ |_  yht  ^ 

We  then  show  that  this  quadratic  form  can  be  bounded  uniformly  from  above  and  below 
by  the  quadratic  form  of  the  composite  finite  element  problem.  Finally  we  establish  that 
we  can  work  practically  with  this  preconditioner  and  that  the  error  propagation  matrix  is 
the  same  as  that  of  the  symmetrized  multiplicative  algorithm,  introduced  in  the  previous 
section.  The  algorithms  are  thus  the  same. 

In  addition  to  the  projections  P/,  we  also  use  another  family  of  operators,  ifj,  i  = 
j  -  l,j,  defined  by 

Hivhix)  =  vh{x),  xen\ Qj, 

a{H}vH,  Wh)  =  0,    \/wh  e  V^'  n  ^o'(".)- 

We  can  call  HjVh  the  /i^  — harmonic  extension  of  v^  to  Qj,  since  it  is  the  solution  in 
V'^'  of  a  discrete  Dirichlet  problem  with  zero  right  hand  side  and  with  boiindary  data  on 
dQj  given  by  Vh- 

In  order  to  prepare  for  the  work  that  remains,  we  formulate  three  lemmas. 

Lemma  1.  H\zIH\-^  =  H\zl,  3  <  i  <  k. 

The  proof  follows  immediately  from  the  definition. 

Lemma  2.  a{Hl-^Uh,Hl~^Vh)  =  a(i7;-^Uh,  i;^),  Vu^  G  V'',  Vi'i,  G  F''-',  2<i<  k. 

Proof:  From  the  definition  oi-H]'^  foUows  that  HI'^vh  -  i-h  G  V^-'  r\H^{Qi).  Since 
Hl~  Uh  is  /ii_i —harmonic  on  fi,  the  result  follows  by  the  orthogonality  inherent  in  the 
definition  of  /ii_i— harmonic  functions. 

The  following  lemma  is  a  consequence  of  the  extension  theorem  given  in  Widlund  [15]. 

Lemma  3.  There  exists  a  constant  C  which  is  independent  of  u^  and  the  mesh  sizes, 
such  that  for  i  =  j  —  1,_/, 

an,{H'jUh,H]uh)  <  CaQ^_^\Q^{uh,Uh),  Vu^  G  V^. 

We  wiU  not  discuss  the  proof  of  this  lemma.  We  note  that  the  constant  C  necessarily 
blows  up  if  we  let  the  area  of  Qj-i  \  Qj  shrink  to  zero,  keeping  Qj-i  fixed.  Such  situations 
are  not  of  particular  interest  in  our  applications. 

Trivially,  we  can  write  Uh  =  Pj^Uh  +  {I  -  Pk)uh-  It  is  easy  to  see  that  the  second  term 
equals  H^Uh.  The  two  terms  are  orthogonal  in  the  sense  of  the  bilinear  form  and  thus 

a{uf,,VH)  =  aiPl:uh,Pl:vH)  +  a{(I  -  P^)uh,{I  -  Pk)vh) 
=  a(P!:uH,Pl:vH)  +  a{H^uu,Htvh)  • 

We  now  introduce  a  preconditioner,  which  in  the  case  k=2  is  the  final  one  but  wliich  in  the 
general  case  is  only  a  first  step  of  our  construction. 

a{Pl:uH,Pl:vH)  +  a{H^-'uf,,H'^''v^)  . 


The  /ifc -harmonic  function  on  fife  has  thus  been  replaced  by  the  /ifc_i -harmonic  function 
with  the  same  boundary  values.  Since  the  latter  space  of  discrete  harmonic  fimctions  is 
smaller,  it  is  easy  to  see  that  the  preconditioner  is  bounded  from  below  by  the  original 
quadratic  form 

aiuh,Uh)  <  a{Pl:uH,P!:uh)  +  aiH^-\H,H^-\h)  ■ 

To  find  a  bound  from  above,  we  have  to  show  that  the  energy  attributable  to  the  hk-i- 
harmonic  function  ccin  be  bounded  by  that  of  the  /ifc -harmonic  function.  This  follows 
directly  from  Lemma  3. 

When  we  now  tiirn  to  the  case  of  fc  >  2,  we  repeatedly  replace  certadn  discrete  harmonic 
functions  by  others.  In  each  step,  we  replace  the  current  last  term  of  the  preconditioner  by 
two.  We  begin  by  considering  the  identity 

It  is  easy  to  see  that  the  second  term  equals  H^zIH^'^uh-  As  in  the  first  step,  we  replace 
this  last  term  by  another,  namely  HJ^ZiH^'^Uh,  and  we  obtain  a  preconditioner  with  three 
terms: 

a{Pi:uH,pi:vh)  +  a{Pi::iHi^-'uH,p;::iHt'vH)+ 

We  simplify  this  expression,  by  using  Lemma  1  and  replace  the  third  term  of  the  precon- 
ditioner by  a{H^Zluh, H^Zl^h)-  By  repeating  the  process  just  outlined,  we  arrive  at  the 
following  preconditioner: 

b{un,VH)  =  a{PJ^UH,Ptvh)  +  a{PlzlHl-'uH.P'^ZlHl-'vu)  ^^^ 

^-■■■a{PlHluH,PlHlvh)  +  a{H\uh,Hlvh). 

We  can  now  complete  our  proof  of  the  optimality  of  the  preconditioner.  By  using  the 
definitions  of  the  operators  P,*  bjxA  Hj,  we  find  that 

a{Pl^Uh,Ptuh)  =  an^{Pl:uh,Pl:uh) 

and,  for  3  <  i  <  fc, 

aiPizlHi-'uH,PiZ^Hi-'uH)  =  ao,_APi:lHi-'uH,PizlHi-'uH) 

=  ao^^ASfuH^Hi-'uH)  -  an,_AHizluh,HizluH). 

In  the  last  step,  we  have  used  Lemma  1. 

By  using  the  definition  of  the  Hj  operators,  we  can  rewrite  the  preconditioner  as 

+  ani_x\ni("h,"/.)  +  an^(H!^~^Uh,H^~^Uh)  -  an^^^{Hf.Zluh,H^Zly-h) 
+  ani_3\ni-,(wh>"h)  +  an^-AHkZi^h,Hk-i^h)  -  an^_^{H^Zluh,H^zluh) 
+  •••ani\n3("h»"h)  +  an^iHluh^HzUh) 
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This  quadratic  form  can  be  written  as  the  sum  of  an(uh,  Wh)  and  a  number  of  terms  of  the 
form 

Since,  on  fi;,  Hl~  Uh  and  H^Uh  are  /ii_i— harmonic  and  /i,  — harmonic  functions,  respec- 
tively, with  the  same  boundary  values  on  dQi,  it  is  easy  to  see  that  these  terms  are  positive. 
This  shows  that  the  preconditioner  is  bounded  from  below  by  aQ{uh,Uh).  To  get  an  upper 
bound,  we  only  have  to  estimate  the  positive  terms.  By  Lemma  3,  the  energy  attributable 
to  the  subregion  Qi,  of  the  /i,_i  —harmonic  function  Hl~  Uh  can  be  estimated  by  a  constant 
times  an,_i\ni("/i)  ^h)-  The  proof  of  the  upper  bound  is  completed  by  summing  over  i. 

What  remains  is  to  establish  that  the  use  of  this  preconditioner  leads  to  a  series 
of  problems  which  are  directly  related  to  the  projections  PI  in  particular  the  symmetric 
multiplicative  algorithm  given  in  section  2.  For  the  unaccelerated  algorithm,  the  equation 
that  we  have  to  solve  is  of  the  form 

b{Suh,Vh)  =  a(eh,Vh)  , 

where  eh  is  the  error  before  the  current  step  and  Suh  the  next  correction.  We  first  choose 
test  functions  Vh  in  the  subspace  V'^^ .  All  terms  of  the  preconditioner,  except  the  iirst, 
vanish  for  such  test  functions.  Therefore  the  result  of  this  first  fractional  step  equals 
P^Suh  =  -Pfc  £;,.  We  choose  F'*'-'  as  our  second  space  of  test  functions.  All  the  terms 
of  the  preconditioner  except  the  first  two  vanish  and  the  first  is  already  known  and  can 
therefore  be  moved  over  to  the  right  hand  side.  A  straight  forward  calculation  and  Lemma 
2  show  that  the  new  right  hand  side  is  equal  to  a((l  -  Pj^)eh,  Vh)-  In  this  second  fractional 
step,  we  obtain  P^Z^  H^~^8uh  =  P^Iiil  -  Pk)<^h-  Proceeding  in  this  manner,  we  compute, 
in  the  k-th.  fractional  step,  H\6uh  =  P/  {I  -  P^)--  ■  (/  -  P^)eh.  At  this  time,  the  values  of 
6uh  are  available  on  fii  \fi2-  We  can  use  its  boundary  values  on  5^2  as  Dirichlet  values  and 
solve  a  problem  in  V^^  to  obtain  the  values  of  8uh  in  ^2  \  ^3.  Proceeding  in  this  manner, 
using  a  total  of  (2fc  -  1)  fractional  steps,  we  finally  obtain  buh  everywhere.  A  calculation 
shows  that  the  error  propagation  operator  corresponding  to  the  whole  step  is  equal  to 

k-\\  IT        T>\\IT         d2\  ,t        nk. 


(/-P,^)(7-P^-/)...(/-P/)(7-P^)...(7-P, 


k 


This  shows  that  the  method  defined  by  the  preconditioner  indeed  is  the  same  as  that  of 
the  symmetrized  multiplicative  algorithm  defined  in  section  2. 

4.  Some  Additive  Algorithms.  We  recall,  that  in  the  so  called  additive  algoritlmis. 
we  solve  equation  (3)  by  the  conjugate  gradient  or  some  other  standard  iterative  algorithm. 
The  different  algorithms  can  simply  be  defined  by  specifying  the  subspaces  V,  or  alterna- 
tively the  projections  P,.  Before  we  discuss  specific  algorithms,  we  make  some  remarks  on 
estimating  the  eigenvalues  of  P  from  above  and  below. 

The  upper  bound  on  the  spectrum  is  obtained  by  bounding 

a{Pvh,Vh)  =  a{PiVh,Vh)  +  a(P2Vh,Vh)  +  •••  +  a[PkVh,Vh)  , 

from  above  in  terms  of  a{vh,Vh).  We  can  use  Schwarz'  inequality  and  the  fact  that  P^  is 
a  projection  to  prove  that  each  term  is  bounded  by  a(vh,Vh)  and  thus  the  spectrmn  of  P 
is  bounded  from  above  by  k.  Our  goal,  however,  is  to  establish  a  uniform  boimd  on  the 
condition  number.  This  can  be  done  if  the  terms  are  orthogonal  or  almost  orthogonal;  cf. 
discussion  below. 


A  lemma  from  Lions  [9]  provides  a  method  for  obtaining  lower  bounds.  Since  the  proof 
of  his  result  is  quite  short,  we  include  it  in  this  paper. 

Lemma  4.  Let  uh  =  Y^J'^i  Uh,i,  where  Uh,i  £  Vi,  he  a  partition  of  an  element  ofV^  and 
assume  further  that  Y,i=i  a(uh,i,Uh,i)  <  C^a{uh,Uh),  Vu;,  e  V".   Then  A„i„(P)  >  Cq^. 

Proof:  By  elementary  properties  of  symmetric  projections  and  the  representation  of 
Uh  as  a  svmi,  we  find  that 


TC  ^  "• 

a{uh,Uh)  =  '^a{uh,Uh,i)  =  ^a{uh,  PiUh,,)  =  2J '^(^i^'i:  ^'H,i) 

Therefore, 


i=i  1=1  1=1 


a{uh,Uh)  <  (^a(P.7i^,P,un))'/'(J^a(t^h,MWh,i))'' 


1/2 

i=l  i=l 

By  the  assumption  of  the  lemma 


a{uh,Uh)  <  ClY,<Pi^h,P^Uh)  =  Cl^a(PiUh,Uh)  =  C^a{Puh,Uh), 

t=l  i=l 

and  the  lemma  is  established. 

We  consider  three  different  algorithms  and  distinguish  between  them  by  using  a  super- 
script 1,  2  or  3.  The  most  natural  algorithm  amoimts  to  using  the  projections  P,  =  P/, 
where  the  projection  operators  P-  have  been  defined  in  equation  (2).  The  condition  number 
of  this  algorithm  grows  Linearly  with  k.  We  have  already  remarked  that  the  eigenvalues  of  P 
always  are  bounded  from  above  by  k.  This  bound  is  attained  if  F'"  nifo(^fc)  i^  not  empty, 
i.e.  when  the  mesh  size  is  fine  enough.  Any  such  fimction  belongs  to  T^  ',  t  =  1,2, . . .,  A;, 
and  is  exactly  reproduced  by  each  of  the  projection  operators.  It  is  thus  an  eigenfunction 
with  eigenvalue  k.  Similarly,  any  function  which  belongs  to  V^'-  H  Hq(Qi  \  ^2)  is  an  eigen- 
fimction  with  eigenvalue  1.  We  will  show  later  that  all  the  eigenvalues  are  boimded  from 
below  by  a  constant,  which  completes  the  proof  that  the  condition  number  of  P^^'  is  of 
order  k. 

A  very  promising  method,  for  which  to  our  knowledge  the  optimality  has  only  been 
established  in  a  quite  special  model  case,  cf.  [10],  is  based  on  using  Pj"  =  P-  -  P/+i  , 
i  <  k  -  1  and  P|.^*  =  P^^  as  the  basic  projections  in  equation  (3).  It  easy  to  show  that  these 
differences  of  projections  are  projections  and  that  the  composite  finite  element  space  V 
is  the  direct  sum  of  the  corresponding  subspaces.  A  difficulty  experienced  when  trying  to 
establish  that  the  eigenvalues  of  P'^*  are  bounded  uniformly  from  below,  by  using  Lions' 
lemma,  is  related  to  this  complete  lack  of  flexibility  in  representing  a  given  element  of 
V^  as  a  simi  of  elements  of  the  subspaces.  So  far  we  have  only  been  able  to  prove  that 
C'o  <  const,  t. 

We  can,  however,  prove  the  following  result. 

Theorem  2.    The  eigenvalues  of  P''"^  are  uniformly  bounded  from  above  by  a  constant. 

Before  we  turn  to  the  proof  proper,  we  observe  that  the  unbalance  of  the  first  method, 
which  resulted  in  the  amplification  of  certain  functions  by  a  factor  k,  is  no  longer  possible. 
By  regrouping  terms,  we  also  see  that 

p(2)  =  p;  +  {p^  -  P2')  +  •••  +  (^fc  -  Pk'')  ■  (5) 

10 


All  the  terms,  except  the  first,  are  similar  to  miiltigrid  corrections  since  they  each  represent 
the  difference  between  two  solutions  on  the  same  subregion,  using  two  different  mesh  sizes. 
By  using  the  representation  of  P*"',  given  in  equation  (5),  we  obtain 

P^'^UH  =  P/  UH  +  (P|   -P^)UH  +  ---+  {Pt  -  P^'')UH 

=  Ui   +  U2   + \-  Uk- 

We  show  that  these  terms  are  increasingly  orthogonal: 

a{u(,Um)  <  const,  ^f    ™'(a(u<,  Ui-))^''^(a(u„,  u„))^/^ 

where  qi  <  I.  uniformly.  This  is  a  so-called  strengthened  Cauchy  Inequality.  Without 
loss  of  generality,  we  assume  that  i  <  m.  It  is  easy  to  show  that  aium,  v^)  =  0,  Vu^  G 
yh„-i  nifgl^m).  i-e-  Um  is  /im -i -harmonic  on  Qrn-  Let  U(  =  Va,c()-'  where  the  ())■'  are 
basis  functions  in  V'^' .  Any  such  basis  function  is  orthogonal  to  Um  if  its  support  does  not 
intersect  di^m  since  then  either  its  support  belongs  to  f2^  \  Qrn  ot  <p^'  G  F''"-'  n  ^o(ri^). 
Thus  the  only  contributions  to  a{um-,  "/)  originate  from  the  triangles  of  F  '  which  intersect 
5f2rn-  Consider  one  such  triangle  T.  The  restriction  of  u<  to  T  is  a  linear  function  and  thus 
has  a  constant  gradient.  It  can  also  be  written  as  a  linear  combination  of  basis  functions 
in  V'''"-'.  The  support  of  most  of  these  do  not  intersect  dflrn  and  therefore 

aT{um,ue)  =  aTns{um,y-i)  <  {aT{Um,Um))^^'{aTns(U(,Ui)f'' 
<  Const.  q["^~'^aT{um,yn)^^^aT{ui,U(y^^  . 

Here  5  denotes  the  union  of  the  smaU  triangles  that  are  next  to  (9ilm.  We  also  use  our 
assumptions  on  the  triangulations,  given  in  section  2,  and  the  fact  that  the  gradient  of 
U(  is  constant  on  T  and  that  therefore  its  contribution  to  the  quadratic  form  is  directly 
proportional  to  the  area  of  integration.  The  proof  of  the  strengthened  Cauchy  inequality 
is  completed  by  summing  over  the  triangles  of  V^' . 
We  now  note  that 


Here 


a(P'^'u,,P<  =  'tz,) 

=    > 

a(ui,Uj) 

<  U^AU. 

/     ' 

q\ 

^l        ■ 

■•  'C\ 

A  =  const.        ^^ 

1 

<1\ 

■  ■  <h  " 

V^f-^ 

k-2 

fc-3 

1x 

;;  1  / 

and 

It  is  easy  to  show  that  the  Euclidean  norm  of  A  is  uniformly  bounded.   Thus 

k  k 

a{P'^^u,P^^^u)  <  \A\i^Yla{ui,u,)  <  const.   ^a(tf„«,)- 

i=l  1=1 

To  complete  the  proof,  we  note  that 

a{ui,Ui)  =  a{PluH  -  Pl^'uH.PluH  -  P'^' u^) 
=  a(PluH-Pl-'ut,,PiuH), 
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since  P/u^  -  P/   ^Uh  is  hi-i  -harmonic.  By  using  elementaxy  properties  of  the  projections, 
we  find  that  a{ui,Ui)  =  a(u,-,u;,).  Therefore 

k  k 

t=i  t=i 

Thus 

a{P^^^Uh,P''^^Uh)  <  const.  a{P^^^Uh,Uh), 

from  which  the  upper  bound  on  P*^^  follows. 

If  we  modify  our  projections,  we  arrive  at  an  algorithm  for  which  we  have  a  proof  of 
optimality. 

(3) 

Theorem  3.   The  additive  method  defined  by  the  projections  P,-      =  PI  —  PIj^^i  ^  ^  k  —  2, 

PI._i  —  Pkli    '^''^d  Pk  Pk   ^'^^  ^  condition  number  which  is  independent  of  k  and  the 

number  of  degrees  of  freedom. 

A  uniform  upper  bound  is  obtained  by  using  Theorem  2  and  the  following  formula. 


p(3)  ^  pi  _  pi  ^p2  _p2  _^,,.^  pk-i  +  p^ 


=  Pl  -Pi+P^  -Pi  +  ---  +  Pt 
+  Pl-Pl+Pl-Pl^...  +  P'^Zl 

_  p(2)  +p(2) 

By  Theorem  2,  there  is  an  upper  bound  for  P^^'  and  the  same  result  also  provides  a  uniform 
upper  bound  for  P^^\  which  is  an  operator  corresponding  to  a  composite  finite  element 
problem  using  k  —  1  levels. 

From  this  argimient  and  the  fact  that  P*^'  is  positive  definite,  it  foUows  that  the 
condition  number  of  P^^^  caumot  be  much  better  than  that  of  P*^'.  A  lower  bound  of  the 
spectrum  of  P^"'*  is  given  by  a  partitioning  formula  and  Lions'  Lemma.  Before  we  give  the 
details,  we  note  that  this  argument  also  provides  a  lower  bound  for  the  spectrum  of  P*^' 
since  the  subspaces  related  to  the  projections  of  that  method  include  those  of  P'^'. 

We  partition  an  arbitrary  Uh  G  V^  as 

k 

( 3)    .  (3 ) 

where  V^      is  the  range  of  P^     .  The  formulas  for  the  u/,  ^  are  given  by 

{Uh  on       Qi  \  fij 

hi  -  harmonic     on       ^2  \  ^3 
0  on      ^3 

and  for  2  <  z  <  fc  —  2,  and 

Uh  -  Uh,i-i  on      Qi  \  Cli+i 

Uh.i  =  {  hi  —  harmonic     on      f2i+i  \  ^^+2 
0  on      fit+2 
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Since  the  construction  of  u^  ;,_i  and  Uh,h  is  straightforward,  we  do  not  provide  details.  We 
note  that  on  each  region  Qi  \  fli+i  only  two  of  the  terms  differ  from  zero  cind  that  Uh,i  =  Uh 
on  dVli+i.  By  applying  the  projection  P-^'  to  Uh,i,  we  find  that  u^,,  6  V-  .  The  proof  is 
completed  by  showing  that 

a(uh,t,"h,!)  <  const.  aa^\u^_^^{uh,Uh). 

This  follows  from  a  variant  of  Lemma  3  and  elementary  considerations.  The  proof  of  the 
lower  bound  is  completed  by  simiming  over  i. 
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